home *** CD-ROM | disk | FTP | other *** search
/ Software Vault: The Diamond Collection / The Diamond Collection (Software Vault)(Digital Impact).ISO / cdr44 / newmat08.zip / NEWMAT7.CPP < prev    next >
C/C++ Source or Header  |  1995-01-16  |  24KB  |  746 lines
  1. //$$ newmat7.cpp     Invert, solve, binary operations
  2.  
  3. // Copyright (C) 1991,2,3,4: R B Davies
  4.  
  5. #include "include.h"
  6.  
  7. #include "newmat.h"
  8. #include "newmatrc.h"
  9.  
  10. //#define REPORT { static ExeCounter ExeCount(__LINE__,7); ++ExeCount; }
  11.  
  12. #define REPORT {}
  13.  
  14.  
  15. /***************************** solve routines ******************************/
  16.  
  17. GeneralMatrix* GeneralMatrix::MakeSolver()
  18. {
  19.    REPORT
  20.    GeneralMatrix* gm = new CroutMatrix(*this);
  21.    MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
  22. }
  23.  
  24. GeneralMatrix* Matrix::MakeSolver()
  25. {
  26.    REPORT
  27.    GeneralMatrix* gm = new CroutMatrix(*this);
  28.    MatrixErrorNoSpace(gm); gm->ReleaseAndDelete(); return gm;
  29. }
  30.  
  31. void CroutMatrix::Solver(MatrixRowCol& mcout, const MatrixRowCol& mcin)
  32. {
  33.    REPORT
  34.    Real* el = mcin.store; int i = mcin.skip;
  35.    while (i--) *el++ = 0.0;
  36.    el += mcin.storage; i = nrows - mcin.skip - mcin.storage;
  37.    while (i--) *el++ = 0.0;
  38.    lubksb(mcin.store, mcout.skip);
  39. }
  40.  
  41.  
  42. // Do we need check for entirely zero output?
  43.  
  44. void UpperTriangularMatrix::Solver(MatrixRowCol& mcout,
  45.    const MatrixRowCol& mcin)
  46. {
  47.    REPORT
  48.    Real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
  49.    while (i-- > 0) *elx++ = 0.0;
  50.    int nr = mcin.skip+mcin.storage; elx = mcin.store+nr; Real* el = elx;
  51.    int j = mcout.skip+mcout.storage-nr; int nc = ncols-nr; i = nr-mcout.skip;
  52.    while (j-- > 0) *elx++ = 0.0;
  53.    Real* Ael = store + (nr*(2*ncols-nr+1))/2; j = 0;
  54.    while (i-- > 0)
  55.    {
  56.       elx = el; Real sum = 0.0; int jx = j++; Ael -= nc;
  57.       while (jx--) sum += *(--Ael) * *(--elx);
  58.       elx--; *elx = (*elx - sum) / *(--Ael);
  59.    }
  60. }
  61.  
  62. void LowerTriangularMatrix::Solver(MatrixRowCol& mcout,
  63.    const MatrixRowCol& mcin)
  64. {
  65.    REPORT
  66.    Real* elx = mcin.store+mcout.skip; int i = mcin.skip-mcout.skip;
  67.    while (i-- > 0) *elx++ = 0.0;
  68.    int nc = mcin.skip; i = nc+mcin.storage; elx = mcin.store+i;
  69.    int nr = mcout.skip+mcout.storage; int j = nr-i; i = nr-nc;
  70.    while (j-- > 0) *elx++ = 0.0;
  71.    Real* el = mcin.store+nc; Real* Ael = store + (nc*(nc+1))/2; j = 0;
  72.    while (i-- > 0)
  73.    {
  74.       elx = el; Real sum = 0.0; int jx = j++; Ael += nc;
  75.       while (jx--) sum += *Ael++ * *elx++;
  76.       *elx = (*elx - sum) / *Ael++;
  77.    }
  78. }
  79.  
  80. /******************* carry out binary operations *************************/
  81.  
  82. static GeneralMatrix*
  83.    GeneralAdd(GeneralMatrix*,GeneralMatrix*,AddedMatrix*,MatrixType);
  84. static GeneralMatrix*
  85.    GeneralSub(GeneralMatrix*,GeneralMatrix*,SubtractedMatrix*,MatrixType);
  86. static GeneralMatrix*
  87.    GeneralMult(GeneralMatrix*,GeneralMatrix*,MultipliedMatrix*,MatrixType);
  88. static GeneralMatrix*
  89.    GeneralSolv(GeneralMatrix*,GeneralMatrix*,BaseMatrix*,MatrixType);
  90. static GeneralMatrix*
  91.    GeneralSP(GeneralMatrix*,GeneralMatrix*,SPMatrix*,MatrixType);
  92.  
  93. GeneralMatrix* AddedMatrix::Evaluate(MatrixType mt)
  94. {
  95.    REPORT
  96.    gm1=((BaseMatrix*&)bm1)->Evaluate();
  97.    gm2=((BaseMatrix*&)bm2)->Evaluate();
  98. #ifdef TEMPS_DESTROYED_QUICKLY
  99.    GeneralMatrix* gmx;
  100.    Try { gmx = GeneralAdd(gm1,gm2,this,mt); }
  101.    CatchAll { delete this; Bounce; }
  102.    delete this; return gmx;
  103. #else
  104.    return GeneralAdd(gm1,gm2,this,mt);
  105. #endif   
  106. }
  107.  
  108. GeneralMatrix* SubtractedMatrix::Evaluate(MatrixType mt)
  109. {
  110.    REPORT
  111.    gm1=((BaseMatrix*&)bm1)->Evaluate();
  112.    gm2=((BaseMatrix*&)bm2)->Evaluate();
  113. #ifdef TEMPS_DESTROYED_QUICKLY
  114.    GeneralMatrix* gmx;
  115.    Try { gmx = GeneralSub(gm1,gm2,this,mt); }
  116.    CatchAll { delete this; Bounce; }
  117.    delete this; return gmx;
  118. #else
  119.    return GeneralSub(gm1,gm2,this,mt);
  120. #endif   
  121. }
  122.  
  123. GeneralMatrix* MultipliedMatrix::Evaluate(MatrixType mt)
  124. {
  125.    REPORT
  126.    gm2 = ((BaseMatrix*&)bm2)->Evaluate();
  127.    gm2 = gm2->Evaluate(gm2->Type().MultRHS());     // no symmetric on RHS
  128.    gm1=((BaseMatrix*&)bm1)->Evaluate();
  129. #ifdef TEMPS_DESTROYED_QUICKLY
  130.    GeneralMatrix* gmx;
  131.    Try { gmx = GeneralMult(gm1, gm2, this, mt); }
  132.    CatchAll { delete this; Bounce; }
  133.    delete this; return gmx;
  134. #else
  135.    return GeneralMult(gm1, gm2, this, mt);
  136. #endif   
  137. }
  138.  
  139. GeneralMatrix* SolvedMatrix::Evaluate(MatrixType mt)
  140. {
  141.    REPORT
  142.    gm1=((BaseMatrix*&)bm1)->Evaluate();
  143.    gm2=((BaseMatrix*&)bm2)->Evaluate();
  144. #ifdef TEMPS_DESTROYED_QUICKLY
  145.    GeneralMatrix* gmx;
  146.    Try { gmx = GeneralSolv(gm1,gm2,this,mt); }
  147.    CatchAll { delete this; Bounce; }
  148.    delete this; return gmx;
  149. #else
  150.    return GeneralSolv(gm1,gm2,this,mt);
  151. #endif   
  152. }
  153.  
  154. GeneralMatrix* SPMatrix::Evaluate(MatrixType mt)
  155. {
  156.    REPORT
  157.    gm1=((BaseMatrix*&)bm1)->Evaluate();
  158.    gm2=((BaseMatrix*&)bm2)->Evaluate();
  159. #ifdef TEMPS_DESTROYED_QUICKLY
  160.    GeneralMatrix* gmx;
  161.    Try { gmx = GeneralSP(gm1,gm2,this,mt); }
  162.    CatchAll { delete this; Bounce; }
  163.    delete this; return gmx;
  164. #else
  165.    return GeneralSP(gm1,gm2,this,mt);
  166. #endif   
  167. }
  168.  
  169. // routines for adding or subtracting matrices of identical storage structure
  170.  
  171. static void Add(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  172. {
  173.    REPORT
  174.    Real* s1=gm1->Store(); Real* s2=gm2->Store();
  175.    Real* s=gm->Store(); int i=gm->Storage() >> 2;
  176.    while (i--)
  177.    {
  178.        *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
  179.        *s++ = *s1++ + *s2++; *s++ = *s1++ + *s2++;
  180.    }
  181.    i=gm->Storage() & 3; while (i--) *s++ = *s1++ + *s2++;
  182. }
  183.    
  184. static void Add(GeneralMatrix* gm, GeneralMatrix* gm2)
  185. {
  186.    REPORT
  187.    Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  188.    while (i--)
  189.    { *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; *s++ += *s2++; }
  190.    i=gm->Storage() & 3; while (i--) *s++ += *s2++;
  191. }
  192.  
  193. static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  194. {
  195.    REPORT
  196.    Real* s1=gm1->Store(); Real* s2=gm2->Store();
  197.    Real* s=gm->Store(); int i=gm->Storage() >> 2;
  198.    while (i--)
  199.    {
  200.        *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
  201.        *s++ = *s1++ - *s2++; *s++ = *s1++ - *s2++;
  202.    }
  203.    i=gm->Storage() & 3; while (i--) *s++ = *s1++ - *s2++;
  204. }
  205.  
  206. static void Subtract(GeneralMatrix* gm, GeneralMatrix* gm2)
  207. {
  208.    REPORT
  209.    Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  210.    while (i--)
  211.    { *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; *s++ -= *s2++; }
  212.    i=gm->Storage() & 3; while (i--) *s++ -= *s2++;
  213. }
  214.  
  215. static void ReverseSubtract(GeneralMatrix* gm, GeneralMatrix* gm2)
  216. {
  217.    REPORT
  218.    Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  219.    while (i--)
  220.    {
  221.       *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
  222.       *s = *s2++ - *s; s++; *s = *s2++ - *s; s++;
  223.    }
  224.    i=gm->Storage() & 3; while (i--) { *s = *s2++ - *s; s++; }
  225. }
  226.  
  227. static void SP(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  228. {
  229.    REPORT
  230.    Real* s1=gm1->Store(); Real* s2=gm2->Store();
  231.    Real* s=gm->Store(); int i=gm->Storage() >> 2;
  232.    while (i--)
  233.    {
  234.        *s++ = *s1++ * *s2++; *s++ = *s1++ * *s2++;
  235.        *s++ = *s1++ * *s2++; *s++ = *s1++ * *s2++;
  236.    }
  237.    i=gm->Storage() & 3; while (i--) *s++ = *s1++ * *s2++;
  238. }
  239.    
  240. static void SP(GeneralMatrix* gm, GeneralMatrix* gm2)
  241. {
  242.    REPORT
  243.    Real* s2=gm2->Store(); Real* s=gm->Store(); int i=gm->Storage() >> 2;
  244.    while (i--)
  245.    { *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; *s++ *= *s2++; }
  246.    i=gm->Storage() & 3; while (i--) *s++ *= *s2++;
  247. }
  248.  
  249. // routines for adding or subtracting matrices of different storage structure
  250.  
  251. static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  252. {
  253.    int nr = gm->Nrows();
  254.    MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  255.    MatrixRow mr(gm, StoreOnExit+DirectPart);
  256.    while (nr--) { mr.Add(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
  257. }
  258.    
  259. static void AddDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  260. // Add into first argument
  261. {
  262.    int nr = gm->Nrows();
  263.    MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
  264.    MatrixRow mr2(gm2, LoadOnEntry);
  265.    while (nr--) { mr.Add(mr2); mr.Next(); mr2.Next(); }
  266. }
  267.  
  268. static void SubtractDS
  269.    (GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  270. {
  271.    int nr = gm->Nrows();
  272.    MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  273.    MatrixRow mr(gm, StoreOnExit+DirectPart);
  274.    while (nr--) { mr.Sub(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
  275. }
  276.  
  277. static void SubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  278. {
  279.    int nr = gm->Nrows();
  280.    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
  281.    MatrixRow mr2(gm2, LoadOnEntry);
  282.    while (nr--) { mr.Sub(mr2); mr.Next(); mr2.Next(); }
  283. }
  284.  
  285. static void ReverseSubtractDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  286. {
  287.    int nr = gm->Nrows();
  288.    MatrixRow mr(gm, LoadOnEntry+StoreOnExit+DirectPart);
  289.    MatrixRow mr2(gm2, LoadOnEntry);
  290.    while (nr--) { mr.RevSub(mr2); mr2.Next(); mr.Next(); }
  291. }
  292.  
  293. static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm1, GeneralMatrix* gm2)
  294. {
  295.    int nr = gm->Nrows();
  296.    MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  297.    MatrixRow mr(gm, StoreOnExit+DirectPart);
  298.    while (nr--) { mr.Multiply(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
  299. }
  300.    
  301. static void SPDS(GeneralMatrix* gm, GeneralMatrix* gm2)
  302. // SP into first argument
  303. {
  304.    int nr = gm->Nrows();
  305.    MatrixRow mr(gm, StoreOnExit+LoadOnEntry+DirectPart);
  306.    MatrixRow mr2(gm2, LoadOnEntry);
  307.    while (nr--) { mr.Multiply(mr2); mr.Next(); mr2.Next(); }
  308. }
  309.  
  310. #ifdef __GNUG__
  311. void AddedMatrix::SelectVersion
  312.    (MatrixType mtx, int& c1, int& c2) const
  313. #else
  314. void AddedMatrix::SelectVersion
  315.    (MatrixType mtx, Boolean& c1, Boolean& c2) const
  316. #endif
  317. // for determining version of add and subtract routines
  318. // will need to modify if further matrix structures are introduced
  319. {
  320.    MatrixBandWidth bw1 = gm1->BandWidth();
  321.    MatrixBandWidth bw2 = gm2->BandWidth();
  322.    MatrixBandWidth bwx(-1); if (mtx.IsBand()) bwx = bw1 + bw2;
  323.    c1 = (mtx == gm1->Type()) && (bwx == bw1);
  324.    c2 = (mtx == gm2->Type()) && (bwx == bw2);
  325. }
  326.  
  327. static GeneralMatrix* GeneralAdd(GeneralMatrix* gm1, GeneralMatrix* gm2,
  328.    AddedMatrix* am, MatrixType mtx)
  329. {
  330.    Tracer tr("GeneralAdd");
  331.    int nr=gm1->Nrows(); int nc=gm1->Ncols();
  332.    if (nr!=gm2->Nrows() || nc!=gm2->Ncols()) 
  333.       Throw(IncompatibleDimensionsException());
  334.    Compare(gm1->Type() + gm2->Type(),mtx);
  335. #ifdef __GNUG__
  336.    int c1,c2; am->SelectVersion(mtx,c1,c2);
  337. #else
  338.    Boolean c1,c2; am->SelectVersion(mtx,c1,c2); // causes problems for g++
  339. #endif
  340.    if (c1 && c2)
  341.    {
  342.       if (gm1->reuse()) { REPORT Add(gm1,gm2); gm2->tDelete(); return gm1; }
  343.       else if (gm2->reuse()) { REPORT Add(gm2,gm1); return gm2; }
  344.       else
  345.       {
  346.          REPORT GeneralMatrix* gmx = gm1->Type().New(nr,nc,am);
  347.          gmx->ReleaseAndDelete(); Add(gmx,gm1,gm2); return gmx;
  348.       }
  349.    }
  350.    else
  351.    {
  352.       if (c1 && gm1->reuse() )               // must have type test first
  353.       { REPORT AddDS(gm1,gm2); gm2->tDelete(); return gm1; }
  354.       else if (c2 && gm2->reuse() )
  355.       { REPORT AddDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2; }
  356.       else
  357.       {
  358.          REPORT
  359.      GeneralMatrix* gmx = mtx.New(nr,nc,am); AddDS(gmx,gm1,gm2);
  360.      if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
  361.          gmx->ReleaseAndDelete(); return gmx;
  362.       }
  363.    }
  364. }
  365.  
  366.  
  367. static GeneralMatrix* GeneralSub(GeneralMatrix* gm1, GeneralMatrix* gm2,
  368.    SubtractedMatrix* sm, MatrixType mtx)
  369. {
  370.    Tracer tr("GeneralSub");
  371.    Compare(gm1->Type() + gm2->Type(),mtx);
  372.    int nr=gm1->Nrows(); int nc=gm1->Ncols();
  373.    if (nr!=gm2->Nrows() || nc!=gm2->Ncols())
  374.       Throw(IncompatibleDimensionsException());
  375. #ifdef __GNUG__
  376.    int c1,c2; sm->SelectVersion(mtx,c1,c2);
  377. #else
  378.    Boolean c1,c2; sm->SelectVersion(mtx,c1,c2); // causes problems for g++
  379. #endif
  380.    if (c1 && c2)
  381.    {
  382.       if (gm1->reuse())
  383.       { REPORT Subtract(gm1,gm2); gm2->tDelete(); return gm1; }
  384.       else if (gm2->reuse()) { REPORT ReverseSubtract(gm2,gm1); return gm2; }
  385.       else
  386.       {
  387.          REPORT
  388.      GeneralMatrix* gmx = gm1->Type().New(nr,nc,sm);
  389.          gmx->ReleaseAndDelete(); Subtract(gmx,gm1,gm2); return gmx;
  390.       }
  391.    }
  392.    else
  393.    {
  394.       if ( c1 && gm1->reuse() )
  395.       { REPORT  SubtractDS(gm1,gm2); gm2->tDelete(); return gm1; }
  396.       else if ( c2 && gm2->reuse() )
  397.       {
  398.          REPORT
  399.          ReverseSubtractDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2;
  400.       }
  401.       else
  402.       {
  403.          REPORT
  404.      GeneralMatrix* gmx = mtx.New(nr,nc,sm); SubtractDS(gmx,gm1,gm2);
  405.      if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
  406.      gmx->ReleaseAndDelete(); return gmx;
  407.       }
  408.    }
  409. }
  410.  
  411. static GeneralMatrix* GeneralMult1(GeneralMatrix* gm1, GeneralMatrix* gm2,
  412.    MultipliedMatrix* mm, MatrixType mtx)
  413. {
  414.    REPORT
  415.    Tracer tr("GeneralMult1");
  416.    int nr=gm1->Nrows(); int nc=gm2->Ncols();
  417.    if (gm1->Ncols() !=gm2->Nrows())
  418.       Throw(IncompatibleDimensionsException());
  419.    GeneralMatrix* gmx = mtx.New(nr,nc,mm);
  420.  
  421.    MatrixCol mcx(gmx, StoreOnExit+DirectPart);
  422.    MatrixCol mc2(gm2, LoadOnEntry);
  423.    while (nc--)
  424.    {
  425.       MatrixRow mr1(gm1, LoadOnEntry, mcx.Skip());
  426.       Real* el = mcx();                              // pointer to an element
  427.       int n = mcx.Storage();
  428.       while (n--) { *(el++) = DotProd(mr1,mc2); mr1.Next(); }
  429.       mc2.Next(); mcx.Next();
  430.    }
  431.    gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
  432. }
  433.  
  434. static GeneralMatrix* GeneralMult2(GeneralMatrix* gm1, GeneralMatrix* gm2,
  435.    MultipliedMatrix* mm, MatrixType mtx)
  436. {
  437.    // version that accesses by row only - not good for thin matrices
  438.    // or column vectors in right hand term. Needs fixing
  439.    REPORT
  440.    Tracer tr("GeneralMult2");
  441.    int nr=gm1->Nrows(); int nc=gm2->Ncols();
  442.    if (gm1->Ncols() !=gm2->Nrows())
  443.       Throw(IncompatibleDimensionsException());
  444.    GeneralMatrix* gmx = mtx.New(nr,nc,mm);
  445.  
  446.    Real* el = gmx->Store(); int n = gmx->Storage();
  447.    while (n--) *el++ = 0.0;
  448.    MatrixRow mrx(gmx, LoadOnEntry+StoreOnExit+DirectPart);
  449.    MatrixRow mr1(gm1, LoadOnEntry);
  450.    while (nr--)
  451.    {
  452.       MatrixRow mr2(gm2, LoadOnEntry, mr1.Skip());
  453.       el = mr1();                              // pointer to an element
  454.       n = mr1.Storage();
  455.       while (n--) { mrx.AddScaled(mr2, *el++); mr2.Next(); }
  456.       mr1.Next(); mrx.Next();
  457.    }
  458.    gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
  459. }
  460.  
  461. static GeneralMatrix* mmMult(GeneralMatrix* gm1, GeneralMatrix* gm2)
  462. {
  463.    // matrix multiplication for type Matrix only
  464.    REPORT
  465.    Tracer tr("MatrixMult");
  466.  
  467.    int nr=gm1->Nrows(); int ncr=gm1->Ncols(); int nc=gm2->Ncols();
  468.    if (ncr != gm2->Nrows()) Throw(IncompatibleDimensionsException());
  469.  
  470.    Matrix* gm = new Matrix(nr,nc); MatrixErrorNoSpace(gm);
  471.  
  472.    Real* s1=gm1->Store(); Real* s2=gm2->Store(); Real* s=gm->Store();
  473.    
  474.    if (ncr)
  475.    {
  476.       while (nr--)
  477.       {
  478.          Real* s2x = s2; int j = ncr;
  479.          Real* sx = s; Real f = *s1++; int k = nc;
  480.          while (k--) *sx++ = f * *s2x++;
  481.          while (--j)
  482.             { sx = s; f = *s1++; k = nc; while (k--) *sx++ += f * *s2x++; }
  483.          s = sx;
  484.       }
  485.    }
  486.    else *gm = 0.0;
  487.  
  488.    gm->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gm;
  489. }
  490.  
  491. static GeneralMatrix* GeneralMult(GeneralMatrix* gm1, GeneralMatrix* gm2,
  492.    MultipliedMatrix* mm, MatrixType mtx)
  493. {
  494.    if ( Rectangular(gm1->Type(), gm2->Type(), mtx)) return mmMult(gm1, gm2);
  495.    else
  496.    {
  497.       Compare(gm1->Type() * gm2->Type(),mtx);
  498.       int nr = gm2->Nrows(); int nc = gm2->Ncols();
  499.       if (nc <= 5 && nr > nc) return GeneralMult1(gm1, gm2, mm, mtx);
  500.       else return GeneralMult2(gm1, gm2, mm, mtx);
  501.    }
  502. }
  503.  
  504. #ifdef __GNUG__
  505. void SPMatrix::SelectVersion
  506.    (MatrixType mtx, int& c1, int& c2) const
  507. #else
  508. void SPMatrix::SelectVersion
  509.    (MatrixType mtx, Boolean& c1, Boolean& c2) const
  510. #endif
  511. // for determining version of SP routines
  512. // will need to modify if further matrix structures are introduced
  513. {
  514.    MatrixBandWidth bw1 = gm1->BandWidth();
  515.    MatrixBandWidth bw2 = gm2->BandWidth();
  516.    MatrixBandWidth bwx(-1);  if (mtx.IsBand()) bwx = bw1.minimum(bw2);
  517.    c1 = (mtx == gm1->Type()) && (bwx == bw1);
  518.    c2 = (mtx == gm2->Type()) && (bwx == bw2);
  519. }
  520.  
  521. static GeneralMatrix* GeneralSP(GeneralMatrix* gm1, GeneralMatrix* gm2,
  522.    SPMatrix* am, MatrixType mtx)
  523. {
  524.    Tracer tr("GeneralSP");
  525.    int nr=gm1->Nrows(); int nc=gm1->Ncols();
  526.    if (nr!=gm2->Nrows() || nc!=gm2->Ncols()) 
  527.       Throw(IncompatibleDimensionsException());
  528.    Compare(gm1->Type().SP(gm2->Type()),mtx);
  529. #ifdef __GNUG__
  530.    int c1,c2; am->SelectVersion(mtx,c1,c2);
  531. #else
  532.    Boolean c1,c2; am->SelectVersion(mtx,c1,c2); // causes problems for g++
  533. #endif
  534.    if (c1 && c2)
  535.    {
  536.       if (gm1->reuse()) { REPORT SP(gm1,gm2); gm2->tDelete(); return gm1; }
  537.       else if (gm2->reuse()) { REPORT SP(gm2,gm1); return gm2; }
  538.       else
  539.       {
  540.          REPORT GeneralMatrix* gmx = gm1->Type().New(nr,nc,am);
  541.          gmx->ReleaseAndDelete(); SP(gmx,gm1,gm2); return gmx;
  542.       }
  543.    }
  544.    else
  545.    {
  546.       if (c1 && gm1->reuse() )               // must have type test first
  547.       { REPORT SPDS(gm1,gm2); gm2->tDelete(); return gm1; }
  548.       else if (c2 && gm2->reuse() )
  549.       { REPORT SPDS(gm2,gm1); if (!c1) gm1->tDelete(); return gm2; }
  550.       else
  551.       {
  552.          REPORT
  553.      GeneralMatrix* gmx = mtx.New(nr,nc,am); SPDS(gmx,gm1,gm2);
  554.      if (!c1) gm1->tDelete(); if (!c2) gm2->tDelete();
  555.          gmx->ReleaseAndDelete(); return gmx;
  556.       }
  557.    }
  558. }
  559.  
  560.  
  561. static GeneralMatrix* GeneralSolv(GeneralMatrix* gm1, GeneralMatrix* gm2,
  562.    BaseMatrix* sm, MatrixType mtx)
  563. {
  564.    REPORT
  565.    Tracer tr("GeneralSolv");
  566.    Compare(gm1->Type().i() * gm2->Type(),mtx);
  567.    int nr=gm1->Nrows(); int nc=gm2->Ncols();
  568.    if (gm1->Ncols() !=gm2->Nrows())
  569.       Throw(IncompatibleDimensionsException());
  570.    GeneralMatrix* gmx = mtx.New(nr,nc,sm); MatrixErrorNoSpace(gmx);
  571.    Real* r = new Real [nr]; MatrixErrorNoSpace(r);
  572. #ifndef ATandT
  573.    MONITOR_REAL_NEW("Make   (GenSolv)",nr,r)
  574.                            // deleted for ATandT, to balance deletion below
  575. #endif
  576.    GeneralMatrix* gms = gm1->MakeSolver();
  577.    Try
  578.    {
  579.       MatrixCol mcx(gmx, r, StoreOnExit+DirectPart);   // copy to and from r
  580.          // this must be inside Try so mcx is destroyed before gmx
  581.       MatrixCol mc2(gm2, r, LoadOnEntry);
  582.       while (nc--) { gms->Solver(mcx, mc2); mcx.Next(); mc2.Next(); }
  583.    }
  584.    CatchAll
  585.    {
  586.       gms->tDelete();
  587.       delete gmx;                   // <--------------------
  588.       gm2->tDelete();
  589. #ifndef ATandT
  590.       MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
  591.                           // ATandT version 2.1 gives an internal error
  592. #endif
  593. #ifdef Version21
  594.       delete [] r;
  595. #else
  596.       delete [nr] r;
  597. #endif
  598.       Bounce;
  599.    }
  600.    gms->tDelete(); gmx->ReleaseAndDelete(); gm2->tDelete();
  601. #ifndef ATandT
  602.    MONITOR_REAL_DELETE("Delete (GenSolv)",nr,r)
  603.                           // ATandT version 2.1 gives an internal error
  604. #endif
  605. #ifdef Version21
  606.    delete [] r;
  607. #else
  608.    delete [nr] r;
  609. #endif
  610.    return gmx;
  611. }
  612.  
  613. GeneralMatrix* InvertedMatrix::Evaluate(MatrixType mtx)
  614. {
  615.    // Matrix Inversion - use solve routines
  616.    Tracer tr("InvertedMatrix::Evaluate");
  617.    REPORT
  618.    gm=((BaseMatrix*&)bm)->Evaluate();
  619.    int n = gm->Nrows(); DiagonalMatrix I(n); I=1.0;
  620. #ifdef TEMPS_DESTROYED_QUICKLY
  621.     GeneralMatrix* gmx;
  622.     Try
  623.     {
  624.         Compare(gm->Type().i(),mtx);
  625.         SkipConversionCheck SCC;              // otherwise inverting
  626.                                                           // symmetric causes a problem
  627.         gmx = GeneralSolv(gm,&I,this,mtx);
  628.     }
  629.    CatchAll { delete this; Bounce; }
  630.    delete this; return gmx;
  631. #else
  632.     Compare(gm->Type().i(),mtx);
  633.     SkipConversionCheck SCC;              // otherwise inverting
  634.                                                       // symmetric causes a problem
  635.    return GeneralSolv(gm,&I,this,mtx);
  636. #endif   
  637. }
  638.  
  639.  
  640. /*************************** norm functions ****************************/
  641.  
  642. Real BaseMatrix::Norm1() const
  643. {
  644.    // maximum of sum of absolute values of a column
  645.    REPORT
  646.    GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  647.    int nc = gm->Ncols(); Real value = 0.0;
  648.    MatrixCol mc(gm, LoadOnEntry);
  649.    while (nc--)
  650.       { Real v = mc.SumAbsoluteValue(); if (value < v) value = v; mc.Next(); }
  651.    gm->tDelete(); return value;
  652. }
  653.  
  654. Real BaseMatrix::NormInfinity() const
  655. {
  656.    // maximum of sum of absolute values of a row
  657.    REPORT
  658.    GeneralMatrix* gm = ((BaseMatrix&)*this).Evaluate();
  659.    int nr = gm->Nrows(); Real value = 0.0;
  660.    MatrixRow mr(gm, LoadOnEntry);
  661.    while (nr--)
  662.       { Real v = mr.SumAbsoluteValue(); if (value < v) value = v; mr.Next(); }
  663.    gm->tDelete(); return value;
  664. }
  665.  
  666. /********************** Concatenation and stacking *************************/
  667.  
  668. GeneralMatrix* ConcatenatedMatrix::Evaluate(MatrixType mtx)
  669. {
  670.    REPORT
  671.    Tracer tr("Concatenate");
  672. #ifdef TEMPS_DESTROYED_QUICKLY
  673.       Try
  674.       {
  675.          gm2 = ((BaseMatrix*&)bm2)->Evaluate();
  676.          gm1 = ((BaseMatrix*&)bm1)->Evaluate();
  677.          Compare(gm1->Type() | gm2->Type(),mtx);
  678.          int nr=gm1->Nrows(); int nc = gm1->Ncols() + gm2->Ncols();
  679.          if (nr != gm2->Nrows()) Throw(IncompatibleDimensionsException());
  680.          GeneralMatrix* gmx = mtx.New(nr,nc,this);
  681.          MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  682.          MatrixRow mr(gmx, StoreOnExit+DirectPart);
  683.          while (nr--) { mr.ConCat(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
  684.          gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); delete this;
  685.          return gmx;
  686.       }
  687.       CatchAll  { delete this; Bounce; }
  688. #ifndef UseExceptions
  689.       return 0;
  690. #endif 
  691. #else
  692.       gm2 = ((BaseMatrix*&)bm2)->Evaluate();
  693.       gm1 = ((BaseMatrix*&)bm1)->Evaluate();
  694.       Compare(gm1->Type() | gm2->Type(),mtx);
  695.       int nr=gm1->Nrows(); int nc = gm1->Ncols() + gm2->Ncols();
  696.       if (nr != gm2->Nrows()) Throw(IncompatibleDimensionsException());
  697.       GeneralMatrix* gmx = mtx.New(nr,nc,this);
  698.       MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  699.       MatrixRow mr(gmx, StoreOnExit+DirectPart);
  700.       while (nr--) { mr.ConCat(mr1,mr2); mr1.Next(); mr2.Next(); mr.Next(); }
  701.       gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
  702. #endif
  703. }
  704.  
  705. GeneralMatrix* StackedMatrix::Evaluate(MatrixType mtx)
  706. {
  707.    REPORT
  708.    Tracer tr("Stack");
  709. #ifdef TEMPS_DESTROYED_QUICKLY
  710.       Try
  711.       {
  712.          gm2 = ((BaseMatrix*&)bm2)->Evaluate();
  713.          gm1 = ((BaseMatrix*&)bm1)->Evaluate();
  714.          Compare(gm1->Type() & gm2->Type(),mtx);
  715.          int nc=gm1->Ncols();
  716.          int nr1 = gm1->Nrows(); int nr2 = gm2->Nrows();
  717.          if (nc != gm2->Ncols()) Throw(IncompatibleDimensionsException());
  718.          GeneralMatrix* gmx = mtx.New(nr1+nr2,nc,this);
  719.          MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  720.          MatrixRow mr(gmx, StoreOnExit+DirectPart);
  721.          while (nr1--) { mr.Copy(mr1); mr1.Next(); mr.Next(); }
  722.          while (nr2--) { mr.Copy(mr2); mr2.Next(); mr.Next(); }
  723.          gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); delete this;
  724.          return gmx;
  725.       }
  726.       CatchAll  { delete this; Bounce; } 
  727. #ifndef UseExceptions
  728.       return 0;
  729. #endif 
  730. #else
  731.       gm2 = ((BaseMatrix*&)bm2)->Evaluate();
  732.       gm1 = ((BaseMatrix*&)bm1)->Evaluate();
  733.       Compare(gm1->Type() & gm2->Type(),mtx);
  734.       int nc=gm1->Ncols();
  735.       int nr1 = gm1->Nrows(); int nr2 = gm2->Nrows();
  736.       if (nc != gm2->Ncols()) Throw(IncompatibleDimensionsException());
  737.       GeneralMatrix* gmx = mtx.New(nr1+nr2,nc,this);
  738.       MatrixRow mr1(gm1, LoadOnEntry); MatrixRow mr2(gm2, LoadOnEntry);
  739.       MatrixRow mr(gmx, StoreOnExit+DirectPart);
  740.       while (nr1--) { mr.Copy(mr1); mr1.Next(); mr.Next(); }
  741.       while (nr2--) { mr.Copy(mr2); mr2.Next(); mr.Next(); }
  742.       gmx->ReleaseAndDelete(); gm1->tDelete(); gm2->tDelete(); return gmx;
  743. #endif
  744. }
  745.  
  746.